import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

# 定义函数
def f(x):
    return x**2 + 1

# 微元法计算面积
def area_by_infinitesimal(a, b, n=1000):
    """
    使用微元法计算曲边梯形面积
    a, b: 积分区间
    n: 划分的小区间数量
    """
    x = np.linspace(a, b, n+1)
    dx = (b - a) / n
    area_elements = f(x[:-1]) * dx  # 每个小微元的面积
    total_area = np.sum(area_elements)
    return total_area

# 积分区间
a, b = 0, 2

# 使用scipy的积分函数验证结果
area_integrate, error = integrate.quad(f, a, b)

# 不同的n值
n_values = [10000, 1000000, 10000000]

for n in n_values:
    # 计算面积
    area = area_by_infinitesimal(a, b, n)
    print(f"n = {n:>10,}:")
    print(f"微元法计算的面积: {area:.6f}")
    print(f"积分函数计算的面积: {area_integrate:.6f}")
    print(f"绝对误差: {abs(area - area_integrate):.6f}")
    print()